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ABSTRACT 

We report on next phase of our study of rotating accretion flows onto black 
holes. We consider hydrodynamical (HD) accretion flows with a spherically sym- 
metric density distribution at the outer boundary but with spherical symmetry 
broken by the introduction of a small, latitude-dependent angular momentum. 
We study accretion flows by means of numerical two-dimensional, axisymmetric, 
HD simulations for variety of the adiabatic index, 7 and the gas temperature at 
infinity, Cqq. Our work is an extension of work done by Proga & Begelman who 
consider models for only 7 = 5/3. Our main result is that the flow properties 
such as the topology of the sonic surface and time behavior strongly depend on 
7 but little on c^o- In particular, for 1 < 7 < 5/3, the mass accretion rate shows 
large amplitude, slow time-variability which is a result of mixing between slow 
and fast rotating gas. This temporal behavior differs significantly from that in 
models with 7^5/3 where the accretion rate is relatively constant and from that 
in models with 7^1 where the accretion exhibits small amplitude quasi-periodic 
oscillations. The key parameter responsible for the differences is the sound speed 
of the accretion flow which in turn determines whether the flow is dominated by 
gas pressure, radiation pressure or rotation. Despite these differences the time- 
averaged mass accretion rate in units of the corresponding Bondi rate is a weak 
function of 7 and Cqo- 

Subject headings: accretion - hydrodynamics - black hole physics - outflows - 
galaxies: active - methods: numerical 



1. Introduction 



Many types of astrophysical objects are powered by black hole (BH) accretion. Radi- 
ation energy produced by accretion can be very high and can explain such dramatic phe- 
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nomena as quasars, powerful radio galaxies, X-ray binaries, and gamma ray bursts (GRBs). 
However, BH accretion does not always result in high radiative output. This is true for both 
stellar BH and super massive black holes (SMBH). In particular, objects with SMBH appear 
to spend statistically most of their time in an inactive phase. Inactive SMBHs are not some- 
thing that one would expect, because these black holes are embedded in the relatively dense 
environments of galactic nuclei. Therefore it is natural to suppose that the gravity due to an 
SMBH will draw in matter at high rates, leading to a high system luminosity. Monitoring of 
X-ray binaries with stellar BHs reveals that these objects often exhibit large time variability 
in the total energy output in the spectral energy distribution. Then one of the main goals 
of any theory of BH accretion is two explain why accretion proceeds through very different 
modes. 

Generally, the radiative output from accretion depends on the mass accretion rate 
and an efficiency factor, rj. In previous papers of this series (Proga & Begelman 2003a,b, 
hereafter PB03a and PB03b, respectively), we studied how physical conditions at large dis- 
tances from SMBH affect Ma in the so-called radiatively inefficient accretion flows (RIAF). 
RIAF with very low 77 and also low Ma have been proposed to explain very low radiative 
luminosities in systems as such Sgr A* (e.g., Ichimaru 1977; Rees et al. 1982; Narayan & Yi 
1994, 1995; Abramowicz et al. 1995; Blandford & Begelman 1999; Sharma et al. 2007 and 
reference therein). In PBOSa, we addressed the issue of how Ma depends on the distribution 
of specific angular momentum, I at large radii assuming that the adiabatic index, 7 = 5/3. 

For high /, corresponding to the circularization radius larger than the last stable orbit, 
gas cannot directly accrete onto a BH, unless some physical mechanism like e.g., viscosity or 
magnetic fields, or both, transports angular momentum towards. For inviscid accretion fiow 
the matter with too high I, either fiows outward due to gas pressure and centrifugal forces or 
accumulates close the central accretor (e.g. Hawley et al. 1984a,b,; Clarke et al. 1985; Chen 
et al 1997; PB03a; Janiuk & Proga 2007). The simulations presented by PB03a illustrate 
a general flow pattern with an inflow in the polar funnel, and an equatorial outflow (Fig. 1 
there). Such flow pattern is induced by the angular momentum distribution simply because 
rotating gas tend to converge toward the equator. If / is high, the gas forms a subsonic very 
dynamic torus with gas flowing out. The mass accretion rate can be signiflcantly smaller 
compared to corresponding the Bondi rate. 

In this work, we address the problem how the properties of the accretion flow presented 
in PB03a, change when different micro-physical properties of the flow are assumed. In 
particular, we explore effects of changing 7 in the polytropic equation of state. In reality, 
accretion flows with different 7's, may correspond to different types of objects or different 
phases of activity. For example, in the weakly active galaxies like e.g. Sgr A*, 7 is usually 
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considered to be around 5/3, because these flows are believed to be radiatively inefficient 
and gas pressure dominated. On the other hand, the large energetic output observed in 
GRBs indicates that an accretion flow must be radiation pressure dominated (e.g., Meszaros 
2006) and a relativistic equation of state is required with 7=4/3. Another example are 
protogalactic disks, which are sometimes considered as being formed by isothermal accretion 
flow (7 ~ 1; e.g.. Mo et al. 1998). Also, the so-called 'high' and 'low' accretion states in 
the X-ray binary systems, may reflect physical proprieties changing in time (for a review see 
Done et al. 2007). 

Our work is a straightforward extension of PB03a's work, who showed results for 7=5/3. 
We extend their models by considering the flows, for 7 ranging between 1 and 5/3. We also 
perform additional calculations with sound speed at infinity c^^oo much smaller than that 
considered by PB03a allowing us to model flows with the Bondi radius as large as that 
estimated in real systems, for example in Sgr A*. Taking advantage of faster computers, 
we perform simulations with higher resolution in the 9 direction, and with much larger 
computational domain in comparison to PB03a. We keep other model parameters such as 
angular momentum distribution as in PB03a. 

For 7 = 5/3, our new simulations are consistent with those presented by PB03a. How- 
ever, we flnd that Mq, the flow structure, dcflncd by the gas density and angular momentum 
distribution, and the sonic surface topology depend on 7. In particular, for 7 = 5/3, Ma is 
nearly constant whereas for 7 = 4/3, Af^ shows stochastic, large-amplitude time- variability. 
For 7 = 1.01, Ma shows small-amplitude periodic changes. In Sec. 2, we describe a general 
set-up of our simulations. In Sec. 3, we present our results. In the last section, we discuss 
the results, and relate them to results found in previous studies. 



2. Method 

2.1. Hydrodynamics 

To calculate the structure and evolution of an accreting flow, we solve the equations of 
hydrodynamics 

§ + .V.v = 0. (1) 
P^ = -Vi' + pV4>, (2) 
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where p is the mass density, P is the gas pressure, v is the velocity, e is the internal en- 
ergy density, and $ is gravitational potential. We adopt an adiabatic equation of state 
P = (7 — l)e, where 7 is an adiabatic index. Our calculations are performed in spheri- 
cal polar coordinates {r,6,(f)). We assume axial symmetry about the rotational axis of the 
accretion flow {9 = 0° and 180°). 

We perform simulations using the pseudo-Newtonian potential $ introduced by Paczyhski 
& Wiita (1980) 

$ = (4) 

r-Rs ^ ' 

This potential approximates general relativistic effects in the inner regions, for a non-rotating 

black hole. In particular, the Paczyhski-Wiita (P-W) potential reproduces the last stable 

circular orbit at r = 3i?5 as well as the marginally bound orbit at r = 2Rs. 

Our standard computational domain is defined to occupy the radial range rj = 1.5 Rs < 
r < To — 1.2 Rb and the angular range 0° < ^ < 180°. We consider models with R'g — 10~^, 
R'g — 10~^, and R'g — 10~^. The r — 9 domain is discretized into zones with 140, 180 and 
220 zones (for R'g = 10~^, R'g = 10~^, and R'g = 10~^ respectively) in the r direction and 
200 zones in the 9 direction. We fix zone size ratios, drk+i/drk — 1.05, and d9i/d9i+i — 1.0 
for 0° < ^ < 180°. 



2.2. Initial conditions and boundary conditions 

For the initial conditions of the fluid variables we follow PB03a and adopt a Bondi 
accretion flow. In particular, we adopt Vr and p computed using the Bernoulli function and 
mass accretion rate for spherically symmetric Bondi accretion with the P-W potential. We 
set Poo = 2.2 X 10^^^ g/cm^ and specify Cg^oo through R'g = Rg/Rs (note that R'g = 2c^^^/c^, 
and Rb = is the Bondi radius). Thus R'g characterizes the gas temperature in our 
simulations. We specify the initial conditions by adopting a non-zero / for the outer boundary 
of the flow. 

We consider a case where the angular momentum at the outer radius depends on the 
polar angle via 

liro,9)=lo{l-\cos9\). (5) 
We express the angular momentum on the equator as 

lo = a/R^-RbcSoo, (6) 

where R'q is the circularization radius on the equator in units oi Rb for the Newtonian 
potential (i.e., GM/r"^ = vVr at r = R'cRb)- 
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The boundary conditions are specified as follows. At the poles, (i.e., 6 = 0° and 180°), 
we apply an axis-of-symmetry boundary condition. At both the inner and outer radial 
boundaries, we apply an outflow boundary condition for all dynamical variables. As in 
PB03a, to represent steady conditions at the outer radial boundary, during the evolution of 
each model we continue to apply the constraints that in the last zone in the radial direction, 
vg = 0, V(f, = l{r, 9)/rsm9, and the density is fixed at the outer boundary at all times. Note 
that we allow Vr to float. 

To solve eqs. ([!])-( [3]) we use the ZEUS-2D code described by Stone & Norman (1992), 
modified to implement the P-W potential. 

3. Results 

Here we present results of ten simulations specified by four different values of 7 (i.e., 
7= 5/3, 4/3, 1.2, & 1.01), and three values of R's (i.e., R'g = 10-^ 10-^ and lO'^). The 
simulations for 7 =1.2 were performed only for R'g = 10~^. Summary of all runs is presented 
in Tab [H The table columns (l)-(8) show respectively the name of the run, the numerical 
radial resolution used in the simulation, the value of R'g parameter, the value of circulariza- 
tion radius in units of Bondi radius R'q, 7 adiabatic index, the end time at which we stopped 
each simulation tf (the time is given in units of the dynamical time at the inner boundary 
tdyn= 595 s at r=1.5 R'g for a mass of a black hole to be M^h = 3.6 x IO^Mq), the maximum 
specific angular momentum 1^°-^ at the inner radial boundary, and a time-averaged value 
of mass accretion rate onto a central object in units of Bondi accretion rate Ma/ Mb at the 
end of the simulation. To keep the same resolution at small radii for all cases, our runs 
with lower R^, were performed with more grid points in radial direction. We followed each 
simulation until the quasi-stationary state was achieved, i.e., when the time averaged mass 
accretion rate and torus properties settled down. 

3.1. Mass accretion rate evolution 

Figs. [H [21 and [3] show the mass accretion rate evolution, for R'g = 10^^, R'g = 10~*, 
and R'g = 10~^, respectively, for different 7. In all cases, after an episode of the spherically 
symmetric inflow for which Ma/ Mb =1, Ma/ Mb decreases and starts to fluctuates around 
some time-averaged level. For 7 = 5/3, the mass accretion rate evolves in a similar way for 
all R'g, i.e., it stabilizes quickly and shows no strong time- variability. However, for 7 = 4/3 
and 1.2, Ma the amplitude of the fluctuations is signiflcant. In particular. Ma can suddenly 
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increase by a factor of ~ 5 and then also suddenly decrease. These flares are a result of 
mixing of high and low angular momentum matter (see next sections). For 7 = 4/3, the 
occasional flares appear for all the values of R'g we explored. We suppose that the flares are 
also typical for flows with 7 = 1.2 for regardless of R'g. 

This strong time variability for intermediate 7, is quite surprising given we explore a 
relatively simple HD case. Another surprising result is that the amplitude and time scale of 
the variability depends on 7. In particular, the flares found in models with 7=4/3 and 1.2 
disappears in models with 7 = 1.01. Instead in these models. Ma exhibits small amplitude 
quasi-periodic oscillations (see next sections), for all R'g values. 

Comparing our results for a fixed 7 but various R'g, we find that the time dependence 
is quite insensitive to R'g. This parameter determines the size of the sonic surface and 
therefore can somewhat affect the time-averaged Ma (see Table 1). 

To summarize our results for Ma, in Fig. IHwe plot the mass accretion rate at the end 
of the simulations as a function of R'g for different 7. The final mass accretion rate is a 
relatively complex function of our two model parameters: for 7 = 5/3, Ma decreases with 
increasing R'g whereas for 7 = 1.01 it increases with increasing R'g. For 7 = 4/3, Ma is 
not a monotonic function of R'g, i.e., for R'g < 10~^ it decreases with increasing R'g but for 
R'g ^ 10~^, it increases with increasing R'g. 

As shown by PB03a, Ma depends on the shape and size of the sonic surface. We interpret 
the complex dependence of Ma on 7 and R'q as a result of complex relation between the size 
of the sonic surface that depends on R'g and 7 and the shape of the surface that depends 
strongly on 7 but does not depend significantly on R'g. Next, we will take a closer look at 
detail properties of the simulations. 

3.2. Effects of 7 and 00 on sonic surface topology and flow properties 

Fig. [5] shows the time evolution of the direction of the poloidal velocity and contours of 
the sonic surfaces (i.e., where the Mach number for the poloidal velocity equals 1: \vp/cs\ = 1) 
for four different 7 (from top to bottom, 7= 5/3, 4/3, 1.2 and 1.01). All panels in Fig. [5] 
correspond to models where R'g = 10~^, and each column in the figure corresponds to the 
same time. 

In this and the remaining part of the paper, we focus on the results obtained for R'g = 
10^^ (runs G, H, I, and J), because of relatively similar qualitatively results obtained for 
i?' < 10-3. 
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For 7 = 5/3 (run J, the top panels in Fig. [5]), the initial spherical sonic surface, corre- 
sponding to the Bondi sonic surface, is very small (r^ ~ 26.8Rs) so that it is hard to see it 
in Fig. O Thus we refer a reader to the top right panel in Fig. 9. When the rotating gas 
reaches the initial sonic surface, the rotation quickly modifies the surface shape: the sphere 
turns into two lobes elongated along the rotational axis ("figure eight" shape). This new 
shape is caused by the slowdown of the gas in the equatorial region due to centrifugal force. 
In fact, the centrifugal force and gas pressure halt the inflow and push the gas outward along 
the equator. The outflowing gas creates a shock with the gas infalling from large radii. The 
shock then propagates outward and eventually leaves the computational domain. One can 
follow the evolution of the shock looking at the panels in Fig. [5] from left to right. In the 
right panel showing the end of the simulation, the shock has already left the computational 
domain. The right panel shows also that the flow lost with time its symmetry with respect 
to the equator (we will return to this point in Section 4). 

For 7 = 4/3 (run I, the second row in Fig. [5]), the sonic surface is larger than for 7 = 5/3 
{ts = 253 Rs instead 26.8 -Rs)- At the beginning of the evolution, the initially spherical sonic 
surface does not evolve much, because rotation at the sonic surface is slow. However, as the 
rotating gas reaches small radii it starts to slow down due to increasing centrifugal force. 
Consequently the second, inner sonic surface forms inside the original one. We refer to 
the initial sonic surface as the outer surface. The inner sonic surface forms because high-/ 
matter cannot accrete directly onto the central BH. In run I, gas is more gravitationally 
bound than in run J and does not flow out but rather accumulates in the inner region 
forming a subsonic turbulent torus. The inner flow shocks the supersonically accreting gas 
before the latter feels its centrifugal barrier. We deflne a torus as an equatorial region with 
high-/ gas and subsonic radial velocity. The torus evolves. Namely its density, gas pressure 
and size increase with time. Subsequently the inner sonic surface grows, too. The two sonic 
surfaces eventually connect as the inner surface reaches the outer sonic surface. When this 
happens, the accumulated matter with high / starts to flow out as in run J. The sonic surface 
topology changes and eventually the sonic surface has the flgure eight shape. We note that 
in run I, the flgure eight shape is achieved in a different way than in run J. The second row 
of panels from the top in Fig. Oshow that the matter flows out, but the shock propagation 
is delayed in comparison run J (the top row of panels). At the end of run I, the flow is much 
more asymmetric than for run J (e.g., in run I, there is a strong outflow toward the 'north' 
direction). 

For 7 = 1.2 (run H, the third row of panels in Fig. the evolution of the sonic surface 
topology is similar to that in run I. The main difference is that the inner sonic surfaces 
evolves on longer time scale in run H than run I. The second left panel of Fig. [5] captures 
a moment soon after the two sonic surfaces merged. A delay in formation of the outflow is 
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caused by the fact that for lower 7, the sound speed is lower and the radius of the initial 
outer sonic surface is larger. At the end of run H, the flow is very asymmetric. 

The simulation for 7 = 1.01 (run G, the bottom row of panel in Fig. shows that again 
two sonic surfaces form as in runs H and I. However, the time scale for two sonic surfaces 
to merge is much longer in comparison to other 7 cases, because a large distant between 
the two surfaces (the outer surfaces is at ~500 Rs) and small sound speed. In run H, we 
do not observe the formation of the outflow and we see no shock propagating outward. The 
time scale of the inner sonic surface growth rate is very long, and we are not able to follow 
it, with 200 points in 6 direction, because of too long simulation time. Although not shown 
here, we performed a test simulation, with 2 times lower resolution in 6 direction. In this 
test run we observe the sonic surface topology change in 7 = 1.01 case too, but after much 
longer time in comparison to other 7s (e.g., the topology in run G changes on a time scale 
150 times longer than in run H!). We see no outflow in the equatorial plane for this set of 
parameters, for very long simulation time, even after the figure eight sonic surface is formed. 

We also performed additional computations for 7=1.01 models with the higher flow 
temperatures corresponding to R'g = 10^^ and 2x 10^^, but keeping the same grid parameters 
as for R'g = 10~^ models. As expected, in these additional runs, the flow evolves faster than 
in run H. In particular, an outflow forms on a reasonable time scale. Generally, we found 
that when we suppress the gas pressure enough by decreasing 7 or decreasing Coo, an outflow 
will take a very long time to develop. A delay in development of an outflow is a stronger 
function of 7 than Cqo- 

Fig. [6] compares the mass flux rates time-averaged at the end of simulations as a function 
of radius for runs I, J, H, and G. The mass flux rate at the given radius r is given by: 

M = r'^ (b pVrdQ (7) 

J An 

where dQ = sm{6)d6d(j), p is the density, and Vr is the radial velocity. To calculate the mass 
inflow rate Mi„ at a given radius, we include only contributions with < 0. Whereas to 
calculate the mass outflow rate Mout at a given radius, we include only contributions with 
> 0. To compute the net mass flow flux Mtot we include all contributions regardless of the 
Vr sign. For run J, I, and H, the Mj„ is larger than Mout-, and Mtot is constant. The latter 
indicates that the three models are in time-averaged steady states at the end of simulations. 
We do not show it here, but we check that the total mass in the computational domain 
does not change in time that again indicates the time-averaged steady state. On the other 
hand, run G does not settles into a steady state as matter continues to accumulate between 
bRs < r < lOORs and Mtot changes with radius. 

Fig. [7] compares the mass flux rates as functions of time at the outer boundary of the 
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grid. We calculate Mj„ and Mout from eq. [7] for a given time, in the same way as shown in 
Fig. [6] but assuming that r=ro. For runs J, I, and H, the mass inflow rate is constant for the 
short time at the beginning of the simulations. The gaps in the early time of simulations are 
caused by the shock transition through the outer boundary. In the later times the Mi^ and 
the Mout oscillate and are anticorrelated. For G model, there is no outflow and Mout and 
Min is constant during the whole simulation time. 

Next we present and describe in more details the structure of the flow for different 
7. Figs. [HI El and [TD] show the 2-D structure of various quantities from the simulations 
for different 7, on different scales (i.e., 20 Rs, 200 Rs, and 1200 Rs, in Fig. El El and [lOl 
respectively), at the end of the simulations. Each of the figures consists of four rows, in 
which we present results for 7=5/3, 4/3, 1.2 and 1.01 (panels from top to bottom, J, I, H, 
and G runs respectively). From left to right the panels show the density temperature maps, 
angular momentum contours, angular velocity contours and direction of the poloidal velocity 
overplotted by the sonic surface shape. 

In all four cases, characterized by different 7 index, a torus forms. We distinguish three 
kinds of tori: gas pressure dominated torus, with 7 = 5/3, 'radiation pressure' dominated 
torus with 7 = 4/3, and for 7 = 1.01, due to weak dependency of pressure on density, the 
torus is rotation dominated. 

Run J here corresponds to run B03fla in PB03a. Our results are in full agreement 
with the previous simulations in all respects. The torus is subsonic, the angular momentum 
distribution in the torus is nearly constant / ~ 0.92 Rsc (see also Fig. [11] described below). 
The matter accretes through the polar funnel, and flows out in the equatorial region. The 
mass accretion rate for 7 ~ 5/3 and R'g = 10~^ is consistent with the previous calculations 
of PB03a with two times lower angular resolution. 

For 7 = 4/3 (run I) the radius averaged / in torus is a little higher than in 7=5/3 
case, and it is around 1.1 of the critical value. The density gradient between polar funnels 
and torus is larger than in run J. The sonic surface is strongly asymmetrical. There is the 
above mentioned outflow in 'north' direction. Note that in Figs. [HI [HI and [TUl the angular 
momentum contours undergo strong compression on the 'north' side of the flow. The contours 
are also compressed in the 'southern' part of the grid. This compression is due to the torus 
movement in 'z' direction which we will discussed in more detail in the next subsection. 

In run H, an equatorial outflow is asymmetric and supersonic at large radii. Fig. [TU| 
illustrates strong mixing of high- and low-/ matter at large distances in the 'southern' part 
of the grid leading to an increase in of the mass accretion rate seen as flares in Figs. [H [2l 
and [3] (thick dashed hues) . 
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For 7 = 1.01 (run G), the inflow is a very symmetric, at large radii. On smaller scales 
(r< 200 Rs), the high-Z turbulent torus (red and yellow region in density map in Fig [9l 
bottom panels) is surrounded by slightly lower density, nearly laminar inflow region (the 
green region in the density maps in the same panel). The high-/ matter flows around and 
supplies the matter to the turbulent torus. The torus expands in the 'z' direction with time 
(see Figs [TO] and [T2l) . Thus to reach the black hole, the low-/ matter has to flow around an 
growing in size an obstacle. As a result, an oblique shock, corresponding to the vq sonic 
surface, forms around torus. The oblique shock surrounds the above mentioned region of 
slightly higher density (green region in density maps). For r < 20 Rs (bottom panels in 
Fig. [8]), the torus shape remains quite flat, and is of high density and angular momentum 
with the latter approaching Keplerian distribution. In this inner region, the low-/ gas flows 
toward the equatorial plane almost parallel to the symmetry axis (note that the angular 
momentum contours outside the torus are parallel to the rotational axis). We check that the 
stream lines of the accreting gas, once they reach the torus they follow the torus shape. 

For r < 5 Rs, the low-/ matter forms an accretion cusp. For some time during at the 
initial evolution, the accretion cusp expand in the 'z' direction, because of the torus cusp 
grows in 6 direction. This growth slightly reduces polar funnel accretion, which dominates 
Ma- However soon, the size of the cusp stabilizes. The effect of cusp expansion can be seen 
in the accretion rate evolutionary curve dip (Fig. [1], thin solid line) around t=2.0 x lO'^. Later 
the torus expansion occurs at larger radii only. 

One of the most intriguing property of all the runs is that the angular velocity inside 
the turbulent torus, has cylindrical distribution even though the radial proflle of rotation 
depends on 7. Fig. [TT] shows radial proflles of / along the equator for models with different 7 
at the end of the simulations. For comparison, Keplerian / distribution and / corresponding 
to constant angular velocity are shown too. 

As found by PB03a, / is constant in the torus for 7 = 5/3. We flnd that / does not 
change much inside the torus also for 7=4/3 and 1.2 though the actual value of / in the torus 
increases with decreasing 7. As several other properties, the rotational proflle for 7 = 1.01 
differs from that for other 7. Namely, in run G, the proflle is almost Keplerian for small 
radii (compare solid and dash-dotted line in Fig. [TT]) . 

A detailed inspection of Fig. [TT] shows that / sharply decreases with radius for very small 
radii in runs H and I . This decrease is due to the fact that the inner flow is asymmetric and 
the gas reaching the BH at the end of these runs is not that of the torus but is of the gas 
with low /. The asymmetry also explains why Imax, listed in Table 1, is higher than / at r; 
in Fig. [m i.e., the maximum accreted / does not have to and is not always at the equator. 
The increase of / with radius at large radii is caused by the outer boundary conditions we 
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use with these simulations (i.e., we reset at Tq to its initial value at each time step). 

3.3. Time-variability and asymmetry of the flow 

We find that the temporal behavior and sonic surface topology of the accretion fiow 
depend on 7. In runs with 7 =4/3, 1.2, and 1.01 two sonic surfaces can exit. The sonic 
surface topology evolves in the similar way in all runs with the three 7's although on different 
time scales (the evolution slows down with decreasing 7). During the early evolution, gas 
with large angular momentum (larger than Icr) is locked in a turbulent torus. The torus 
expands but there is no strong outflow because the gas motion is suppressed by supersonically 
infalling matter. During this phase, an oblique shock (seen as a jump in vq in Figs. [8] and 
[9]) forms around the torus. By contrast, in runs with 7 = 5/3, only one sonic surface exists 
and its shape turns from initially spherical to the figure eight like one. The differences in 
the sonic surface topologies are reflected in the mass accretion rate curve: for runs with one 
sonic surface the curves are smooth whereas for runs with two sonic surfaces the curves are 
not smooth (compare various curves in Figs. [H [21 and [3]). 

Generally, the time behavior of the flow depends mainly on the Mach number which in 
our simulations depends on 7 and Coo. To connect the time dependence of Ma with the flow 
properties we start with an analysis of run G with the highest Mach number of the inner 
flow. In run G, shows oscillations which one can attribute to oscillations of the inner 
flow in the z direction (see Fig. [T2l showing time sequence from runs G and H). 

In run G (bottom panels in flg. the high-/ gas that accretes supersonically in the 
equatorial plane from large radii eventually tries to join the turbulent subsonic torus at small 
radii. The torus acts in this case as an obstacle with which the supersonic gas collides at 
about 120 Rs in the equatorial plane. There a strong shock forms. At smaller radii, the 
post-shock gas splits into two streams that flow around the torus. However, the split is not 
into two exactly equal halves. At the early phase of the evolution asymmetry is very small 
as it is due to random numerical asymmetry. However with time the asymmetry grows due 
to strong shocks. In our simulations strong shocks are handled by artiflcial viscosity (Stone 
& Norman 1982). 

The asymmetry amplifled by shocks propagate into other parts of the flow. In particular, 
the torus becomes asymmetric. Subsequently, a larger fraction of the equatorial inflow can 
flow in one of the hemisphere compared to the other hemisphere (Fig. [121 left bottom panel, 
illustrates a phase where more gas flows in the 'southern' hemisphere). Once above the 
torus, the equatorial inflow pushes the polar inflow away from the torus in the 'z' direction 
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(seen as an oblique shock that oscillates and expands in the z direction). The asymmetric 
interaction of the high-Z equatorial and low-/ polar inflows can be also seen in the angular 
momentum contours which are more compressed on one side of the torus. The green colored 
region in Fig. [T2] expands and the density in the region at radii H AORs on the 'southern' 
side of the torus slightly decreases (Fig. [121 middle bottom panel). When this happens, the 
polar inflow of low-/ matter starts to push the shock back toward the equator. The oblique 
shock swings and more of the equatorial inflow goes around the torus from the 'northern' 
side. The scenario repeats (Fig. [121 right bottom panel). As the oscillations continue, the 
oblique shock expands because matter has to flow around the constantly growing torus. 

For 7 = 1.01, the torus is conflned by the supersonic inflow. Therefore, asymmetry is 
strong only at small radii whereas at large radii the flow is symmetric as it is not affected 
by the presence of the asymmetric obstacle. Thus the low amplitude oscillations in the mass 
accretion rate curve are caused by the oscillating torus. The relevant time scale here is the 
dynamical time scale at the radii where the equatorial inflow is split into two streams. Fig. [T3] 
presents a part of the mass accretion curve for model G in the later times of evolution to show 
the periodic variability. The time period between the peaks, corresponds to the dynamical 
time scale at radius of about 84 Rs that corresponds to the outer radius of the torus where 
the high-/ gas splits up into two parts. 

We also calculated a power spectrum using the mass accretion rate curve in the later 
times of evolution. The power spectrum shows broad peak, at frequency corresponding the 
dynamical time scale at radius of about 84 Rs- The peak is broad because the torus increases 
in size with time, in particular the torus radius increases. We computed additional models 
for 7 = 1.01 with higher resolution in 6 direction (400 and 800 points in 6 angle). We found 
that oscillations appear to be independent of the grid resolution. 

For 7 > 1.01, the flow behaves in a similar ways as that in run G. However, there 
is an additional complexity due to high gas pressure and mixing. For 7 > 1.01, the flow 
is less gravitationally bound than in run G and the torus is not as much conflned by the 
supersonic inflow. Therefore, in runs with higher gas pressure the flow is subsonic in a 
relatively big region where mixing between low- and high-/ occurs. Because of the shock 
amplifled asymmetry the mixing leads to large scale asymmetry (compare run G and H in 
Fig. [HI). 

As we mentioned above, for 7 = 4/3 and 7 = 1.2, we also observe occasional bursts 
in the Ma evolution (see e.g. Fig. [T]). The mass accretion rate can signiflcantly increase, 
because of mixing of the high- and low-/ gas at large radii. The flaring behavior of the mass 
accretion rate in models for 7 = 4/3 and 7 = 1.2 is caused by the torus movement in the 'z' 
direction. Fig. [10] (panels corresponding to 7=1.2 case) captured a situation when the torus 
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moved toward the southern hemisphere and partially blocked an inflow of low-/ gas in the 
polar funnel. Then this low-/ gas is pushed toward equator. When the torus later swings 
toward the northern hemisphere, a channel for the low-/ widens one sees a flare in the mass 
accretion rate curve. The reason we do not observe strong flares for 7 = 1.01 is that there is 
no mixing of gas with high- and low-/ at larger distances. 

We find that the flow oscillates in the 'z' direction and is asymmetric with respect to 
the equatorial plane. Instabilities in the accretion disks were studied by many authors. 
Pulsational instability of quasi-adiabatic and quasi-inviscid accretion disks, was examined 
by in many studies (e.g., Kato 1978 Muchotrzeb-Czerny 1986; Kato et al. 1988; Okuda et al. 
1993). Kato et al. (1978) found that the condition for pulsation instability depend on how 
viscosity changes during pulsation. In our simulation only artificial viscosity is introduced, 
which is small. Nevertheless this viscous proves to be important. 

In summary, we find that the asymmetry of the flow is initially induced by the numerical 
effects, but the propagation and amplification of the asymmetry strongly depends on the 
physical conditions in the accretion flow. In the cases where the gas pressure is reduced 
because of low 7 the asymmetry is stronger. For higher 7, the gas pressure is larger and 
the information about the perturbation propagates faster. The gas pressure works against 
effects of shocks and tries to restore the symmetry. We explore effects of gas pressure to 
confirm its role in maintaining symmetry in the flow. To this end, we performed additional 
simulations for the same parameters as run J, but with higher Cqo corresponding to models 
with R'g = 10~^. Fig. [13] shows time evolution of velocity field of one of our exploratory 
runs. The flow with very high Coo is symmetric as one would expect: the information of 
perturbances in high pressure gas propagates vary fast independently of grid effects. This 
simulations supports our reasoning that the asymmetry growing in the flow is not caused by 
numerical effects but rather by shocks amplifying initial seed asymmetry. 

Our time dependent simulations occassionally show structures that resemble features of 
steady state disk solutions, such as backfiows close to the equatorial plane (Kluzniak & Kita, 
2000, and our Fig. [T4l) . or tori in the innermost part of the flow, whose angular momentum 
changes from sub- to super-Keplerian (Figs. [TT], [T^ . Such tori may be relevant to high 
frequency QPOs (e.g. Blaes et al., 2007). 

4. Conclusions 

We report on next phase of our study of rotating accretion flows onto black holes. 
We consider hydrodynamical (HD) accretion flows with a spherically symmetric density 
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distribution at the outer boundary but with spherical symmetry broken by the introduction 
of a small, latitude-dependent angular momentum. We study accretion flows by means of 
numerical two-dimensional, axisymmetric, HD simulations for variety of the adiabatic index, 
7 and the gas temperature at infinity, Cqo- Our work is an extension of work done by PB03a 
who consider models for only 7 = 5/3. We also reran some models from PBOSa with higher 
resolution. The latter are fully consistent with lower resolution PBOSa's runs. 

Our main result is that the flow properties such as the topology of the sonic surface and 
time behavior strongly depend on 7 but little on Cqo- In particular, for 1 < 7 < 5/3, the 
mass accretion rate shows large amplitude, slow time- variability which is a result of mixing 
between slow and fast rotating gas. This temporal behavior differs significantly from that in 
models with 7 ?a 5/3 where the accretion rate is relatively constant and from that in models 
with 7 f=:i 1 where the accretion exhibits small amplitude quasi-periodic oscillations. The key 
parameter responsible for the differences is the sound speed of the accretion flow which in 
turn determines the strength of shocks and whether the flow is dominated but gas pressure 
(7 ^ 5/3), radiation pressure (1 < 7 < 5/3) , or rotation (7 ^ 1). Despite these differences, 
the time-averaged mass accretion rate in units of the corresponding Bondi rate is a weak 
function of 7 and c^. 

We realize that our simulations do not capture several physical processes that may 
be important in real systems (e.g., magnetic flelds, radiative cooling and heating, energy 
dissipation). However, our simulations are done within a general framework so that some 
of these processes can be straightforwardly added (e.g., magnetic fields as in PB03b, or 
radiative processes as in Proga 2007). Therefore, this work could serve as a good reference 
point to analyze and interpret more complete and complex simulations (e.g., we already reran 
some of the models from this paper including magnetic fields, Moscibrodzka and Proga, in 
preparation) . 
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Table 1: Summary of parameter survey. 



Radial 



Run 


resolution 


R'g 






7 


tf 


jmax 
'a 


MJMb 


A 


220 


10-5 


5 X 10- 


-4 


1.01 


1.97 X 10^ 


0.77 


0.45 


B 


220 


10-5 


5 X 10- 


-4 


4/3 


3.10 X 10* 


0.9 


0.25 


C 


220 


10-5 


5 X 10- 


-4 


5/3 


3.14 X 10* 


0.92 


0.13 


D 


180 


10-4 


5 X 10- 


-3 


1.01 


2.80 X 10* 


0.82 


0.40 


E 


180 


10-* 


5 X 10- 


-3 


4/3 


2.75 X 10* 


0.85 


0.13 


F 


180 


10-* 


5 X 10- 


-3 


5/3 


3.19 X 10* 


0.87 


0.15 


G 


140 


10-3 


5 X 10- 


-2 


1.01 


3.48 X 10* 


0.9 


0.30 


H 


140 


10-3 


5 X 10- 


-2 


1.2 


3.10 X 10* 


0.94 


0.20 


I 


140 


10-3 


5 X 10- 


-2 


4/3 


3.56 X 10* 


0.9 


0.17 


J 


140 


10-3 


5 X 10- 


-2 


5/3 


2.67 X 10* 


0.87 


0.30 
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Fig. 1. — Mass accretion rate as a function of time for R'g — 10~^. Time is given in the units 
of dynamical time at the inner boundary. Note the strong flaring behavior for 7 =1.2, and 
4/3. 
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Fig. 2.— As in Fig. 1 but for R'g = 10 
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Fig. 3.— As in Fig. 1 but R's = 10"^ 
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Fig. 4. — Mass accretion rate at the end of the simulations as a function of R'g for different 
7- 
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Fig. 5. — Sequences of poloidal velocity fields, and sonic surface at five different times (from 
left to right t = 177.3, 701.0, 1169.1, 1754.3, and t/ respectively, in units of dynamical time 
at inner radius of computational domain) for 7=5/3, 4/3, 1.2, and 1.01 (from top to bottom). 




Fig. 6. — Mass flux rates as functions of radius for models with R'g — 10~^ and different 
7 time-avcragcd at the end the simulations. For 7=5/3, 4/3, and 1.2, the average is taken 
over A i = 1.75 x 10^ - tf time period. For 7=1.01, we average over A i = 2.5 x 10^ - tf. 
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Fig. 8. — Two-dimensional structure of various quantities at the last time step for 7 ( 
7 index 5/3, 4/3, 1.2, and 1.01 from top to bottom). The panels from left to right are 
snapshots of log p, log T, I, log Q, and of the fluid poloidal velocity direction overplot with 
sonic surface contour. We do not draw labels for angular momentum and angular velocity 
because of strong compression of contours. The angular momentum contours are between 
0.2 and 1.1 with the step of 0.1. The contorus for angular velocity are between -2.5 and -1.5 
with the step of 0.25. The scale of the figure is 20 Rs- 
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Fig. 9. — As in Fig. 8 but on the scale of 200 Rs- The / contours are between 0.2 and 1.1 
with step of 0.1 while the log Q contours are between -5.5 and -3.0 with step of 0.5. 
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Fig. 10. — As in Fig. 8 but on the scale of 1200 Rs (showing the whole computational 
domain). The / contours are between 0.2 and 1.1 with step of 0.1. The log Q contours are 
between -5.5 and -3.0 with step of 0.5. 
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Fig. 11. — Radial profile of angular momentum along the equatorial plane for different 
7 at the end of simulations. For comparison Keplerian angular momentum and angular 
momentum corresponding to constant angular velocity are shown, too. 
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Fig. 12. — Sequence of density maps overplotted by the direction of the poloidal velocity 
for 7=1.2 (upper panels) and 7 = 1.01 (lower panels) for times t=1.97 x 10^, 2.1 x 10^, and 
2.8 X 10^ (upper panels) and for t=9.7 x 10^, 1.12 x 10^, and 1.24 x 10^ (lower panels). These 
time sequences illustrate the evolution of oblique shocks forming in the infalling gas in the 
polar funnel caused by the torus. 
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Fig. 13. — Mass accretion rate curve for model G in the later time of evolution. The curve 
shows a periodic variability with the period corresponding to the dynamical time scale at 
radius of about 84 Rs, where we observe the strongest oscillations the torus and asymmetry 
in the flow (see the lower panels in Fig. 12.) 
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Fig. 14. — Sequence of poloidal velocity fields for five different time for with R'g = 10"^. 
The computational domain and initial condition are the same as in J model, except that the 
Cg^oo is much higher. 



